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Abstract 

Motivated by recent experimental studies, we derive and analyze a two-dimensional model for the contraction patterns 
observed in protoplasmic droplets Physarum polycephalum. The model couples a description of an active poroelastic two- 
phase medium with equations describing the spatiotemporal dynamics of the intracellular free calcium concentration. The 
poroelastic medium is assumed to consist of an active viscoelastic solid representing the cytoskeleton and a viscous fluid 
describing the cytosol. The equations for the poroelastic medium are obtained from continuum force balance and include 
the relevant mechanical fields and an incompressibility condition for the two-phase medium. The reaction-diffusion 
equations for the calcium dynamics in the protoplasm of Physarum are extended by advective transport due to the flow of 
the cytosol generated by mechanical stress. Moreover, we assume that the active tension in the solid cytoskeleton is 
regulated by the calcium concentration in the fluid phase at the same location, which introduces a mechanochemical 
coupling. A linear stability analysis of the homogeneous state without deformation and cytosolic flows exhibits an 
oscillatory Turing instability for a large enough mechanochemical coupling strength. Numerical simulations of the model 
equations reproduce a large variety of wave patterns, including traveling and standing waves, turbulent patterns, rotating 
spirals and antiphase oscillations in line with experimental observations of contraction patterns in the protoplasmic 
droplets. 
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Introduction 

The true slime mold Physarum polycephalum is an extensively 
studied system in biophysics. The plasmodial stage is of particular 
interest, since it exhibits, despite the relatively simple organization 
of this unicellular organism, seemingly "intelligent" physiological 
processes [1]. In this context the term "intelligent" means that, 
given an external stimulus, the plasmodium optimizes its cell 
shape, vein network and growth with respect to transport 
efiiciency as well as robustness with respect to link deletion and 
avoidance of unfavorable conditions [2,3]. Recent experiments 
along these lines show that plasmodia were able to reproduce 
public transport networks on the scale of a petri dish [4] and to 
"solve" maze problems such as finding the shortest path between 
two food sources placed at the exits of a labyrinth [5]. Several 
groups have also investigated the topology and dynamical 
evolution of the vein network in large Physarum plasmodia with 
graph theoretical and statistical physics tools [6-8]. A second 
remarkable phenomenon is the synchronization of the contraction 
patterns in the tubular vein network that generates shuttle 
streaming to distribute nutrients efficiently throughout the 
organism [9] . From the perspective of biophysics it is natural to 
consider these phenomena in the framework of self-organized 



complex systems [10]. For the formulation of mathematical 
models a basic understanding of chemical and mechanical 
processes in the protoplasm is needed. 

A first model for strand contraction combined the viscoelastic 
properties of the ectoplasmic wall with a reaction kinetics that 
regulates the contractile tension of the actomyosin system [11,12]. 
Later, several models in the form of reaction-diffusion [13] and 
reaction-diffusion-advection equations [14,15] were formulated 
that use homogenized quantities, for instance the average strand 
thickness. These models describe Physarum protoplasm as an 
oscillatory medium and treat the mechanical feedback in a 
simplified, qualitative way. More realistic models consider, instead, 
a two-phase description that distinguishes a fluid sol (= cytosol) 
and a solid gel (= cytoskeleton) phase. Some of these models 
account for sol-gel transformations and were used to explain flow- 
channel formation [16] and front dynamics [17]. 

Experiments with microplasmodia, i.e. small plasmodia of sizes 
ranging from 100 /im to several millimeters, provide a possibility to 
study internal amoeboid dynamics of Physarum without the 
pronounced vein structures usually present in Physarum cells of 
larger size. Such microplasmodia are produced by extracting 
cytosol from a Physarum vein and placing it on a substrate. Given 
a sufficient amount of cytosol, protoplasmic droplets will 
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reorganize and form a new independent cellular entity. During the 
first hours of this process such cells show a surprising wealth of 
spatiotemporal mechanical contraction patterns [18,19] (and U. 
Strachauer & M J.B. Hauser, unpublished data). The fact that the 
cell morphology does not change dramatically and that the cell 
does not migrate during the first hours, permits observation of the 
mechanical deformation patterns and waves in a quasi-stationary 
setting. The observed patterns include spirals, traveling and 
standing waves as well as antiphase oscillations (see Fig. 1). 

Various patterns were reproduced previously by a qualitative 
particle-based model [20] . However, this description provided no 
information about the mechanical quantities that are essential to 
understand the intracellular deformation waves and patterns seen 
in the experiments. 

In a more general context, studying the spatiotemporal 
instabilities and the related symmetry breaking in intracellular 
processes has become important to understand many biological 
processes. In a pioneering paper [21], Turing suggested that the 
interplay of reactions and diffusion processes provides a funda- 
mental mechanism for morphogenesis. Later, models for intracel- 
lular pattern formation that include mechanical forces and the 
resulting advection processes have been suggested [22,23]. Active 
gel models describe the cytoskeleton as an active viscous fluid [24] . 
In contrast, experiments on inhomogeneous hydration in cells, 
where large pressure gradients in the cell are observed [25] 
indicate that the cytoplasm can behave like a porous elastic 
sponge-like solid (cytoskeleton) penetrated by a viscous fluid phase 
(cytosol) [26,27]. Moreover, several multiphase flow models have 
been proposed as appropriate descriptions of cytoplasmic dynam- 
ics [28,29]. 

In this paper, we derive and investigate a poroelastic two-phase 
model of the cytoplasm assuming a viscoelastic solid phase 
(cytoskeleton) and a fluid phase (cytosol). Furthermore, we 
incorporate an active tension in the solid phase which is regulated 
by the concentrations of free calcium ions in the fluid phase, that 
are in turn advected by the cytosolic flow field. To account for the 
calcium oscillations observed in experiments with Physarum, a 
simple one-dimensional active poroelastic model derived earlier 
[30] is extended to two dimensions and supplemented by a 
coupling to an oscillatory reaction-diffusion dynamics of the 
intracellular calcium concentration [31]. Unlike the regulating 
species in [30], the total free calcium concentration in the model 




Figure 1. Contractions patterns. Experiments with protoplasmic 
droplets of Physarum polycephalum [18,19]. The color represents the 
local phase of oscillation obtained by a Fourier transformation of the 
spatiotemporal height data: a) standing wave, b) many irregular spirals, 
c) traveling wave, d) antiphase patterns, and e) single spiral. 
doi:1 0.1 371 /journal. pone.0099220.g001 



presented here is not conserved anymore and varies strongly in 
time. As a result, the model displays a short-wavelength instability 
as opposed to the long-wavelength instability found in [30]. The 
choice of the poroelastic approach is motivated by the fact that the 
typical oscillation period of 1 - 2 minutes connected with the 
spatiotemporal contraction patterns discussed above is consider- 
able shorter than the experimentally observed time of 3 - 6 
minutes at which the cytoskeleton starts to exhibit fluid behavior 
[32]. Hence, the resulting model describes the cytoskeleton as an 
active viscoelastic solid coupled to a passive fluid in contrast to 
earlier works that had modeled the cytoskeleton itself as an active 
fluid [22] addressing long time scales, for which fluidization of the 
cytoskeleton has already occurred. 

The inclusion of the calcium oscillator is necessary because it is 
known to be essential in the regulation of the contractile 
actomyosin system of Physarum polycephalum [33]. The mechanical 
and biochemical model parameters are mostly taken from the 
experimental literature on Physarum and therefore allow quanti- 
tative comparisons between model results and experiments. 
Altogether, in this article we derive and analyze a model for the 
intracellular dynamics of protoplasmic droplets that treats the 
cellular mechanics in the framework of a continuous two-phase 
active poroelastic model coupled to an oscillatory biochemical 
medium. Active mechanical stresses induce internal pressure 
gradients that generate cytosolic flow. To account for the latter we 
extend the equations for the internal calcium oscillator proposed in 
[31] to a reaction-diffusion-advection (RDA) model closing the 
mechano-chemical feedback loop. In the methods section we 
introduce and derive the mathematical mode with a description 
divided into a mechanical and a biochemical part. Subsequently, 
the choice of the model parameters is discussed and the numerical 
methods used to discretize and solve the resulting partial 
differential equations (PDEs) are described. The next section 
contains the results obtained by linear stability analysis at the 
homogeneous steady state (HSS) and a two-parameter phase 
diagram with numerical simulations. We present also a selection of 
qualitatively different contraction patterns obtained from simula- 
tions of our model, compare them to earlier experimental findings 
and demonstrate that the variety of patterns found in the 
experiments with Physarum droplets is reproduced successfully. 
The paper is concluded with a discussion of the presented model, 
its limitations and possible extensions. 

Methods 

Model: Mechanical part 

Physarum protoplasm contains a rudimentary form of an 
actomyosin system that is also present in cells of higher vertebrates. 
In contrast to muscle cells the actin filaments in Physarum are 
randomly oriented in the cortex [34,35]. In our model, we assume 
that the cytoplasm contains a solid filamentous phase (gel phase) 
that has viscoelastic properties and exerts contractile tension on 
the system. This is illustrated schematically in Fig. 2. The 
derivation given below is analogous to a recently published 
generic model for active poroelastic media [30]. Unlike the 
regulating species in Ref [30], the total free calcium concentration 
in the model presented here is not conserved anymore and varies 
strongly in time. The fluid part of the cytoplasm is modeled as a 
passive fluid (sol phase) that permeates the cytoskeleton [25]. The 
velocity field in the sol phase will be expressed by the variable v. 
Typical Reynolds numbers that arise from the cytoplasmic flow 
are small [Re«l ) and inertia is negligible. The relative volume 
fraction of solid material is denominated as p^^/ and the fraction of 
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Schematic representation: 




Figure 2. Schematic representation of the the two-phase 
model. Drawing of a Physarum microdroplet (top) in side view 
showing the plasmalemma invaginations filled with extracellular matrix 
(light blue), the fluid phase of the cytoplasm (blue) and the solid 
filamentous phase (black). Top view of the droplet in the simplified 
framework of our two-phase model (bottom). Deformations of the 
cytoskeleton are represented by a distorted grid, flow field in the 
cytosol by arrows and free Ca^^ concentration in yellow, respectively. 
doi:1 0.1 371 /journal. pone.0099220.g002 

the fluid material as p^^i with the additional constraint 

Psol + Pgel = 1 • 

We define a body-reference coordinate system x and a 
displacement field u that gives the deviation of the deformed 
coordinates X: u(x,t) = X(x,t) — x at a time t. The gel velocity is 
the substantial time derivative u = dtU-\-(yu)-x of the displacement 
field. Since the gel is fixed in the reference coordinate system 
(i: = 0), the substantial time derivative u is identical to the partial 
time derivative dfU. 

To determine the flow and displacement field we consider force- 
balance equations of the form 

V-(P.e/KeT + rj))+/gel=0 (1) 



of the medium is random [35]. Hence, we use an isotropic 
constitutive law for the elastic stress-strain relation. It has been 
suggested to consider the cytoplasm as an incompressible medium 
[36]. In the three-dimensional bulk the total mass flux must be 
zero. For small strains this can be expressed as [29]: 

V-(p^e/« + P.o/V) = 0. (3) 

In the following we assume constant sol and gel fractions p^^i 
and p^^i throughout the medium. This is justified, because we 
consider only small deformations. As a result the transport of 
cytosol and the related potential inhomogeneities of the fields p^^/ 
and pggi lead only to second-order corrections in the mechanical 
equations (for details see Text SI in File S2). 

We include a hydrostatic pressure p into the stress tensors of sol 
and gel that originates from the incompressibility of the material 
expressed by Eq. (3) and neglect the osmotic pressures caused by a 
difference in the chemical potential (see Text SI in File S2). We 
assume a Kelvin- Voigt viscoelastic constitutive law (see e.g. [37]) 

k 

for the gel phase. Using Darcy's law w= \p the following 

relation for the drag force is obtained 

/gel = -/sol = pLplelP^^ (4) 

where the parameter is the ratio between the dynamic viscosity rj 
of the cytosol and the permeability k of the porous medium. 
Instead of the fluid velocity v in the laboratory frame, we need to 
consider the velocity w = v — u in the body-reference frame 
introduced above. The sol phase is considered as a passive 
Newtonian fluid. Hence, only viscous stresses are included for the 
sol phase. As a result, Eq. (2) corresponds to the Brinkman 

k 

equation [38]. With the usual expression w= \p for Darcy's 

n 

law and Eq. (2) one can relate the coefficient in our model to the 
viscosity f] of the cytosol and the permeability k of the medium: 
t^p^g^l = r]/k. We divide the overall passive stresses in a dissipative 
and nondissipative (elastic) part: 



^sol- ^sol ^ ^ ^^sol ^ ^ 



V-(p../<oT-)+/soi = 0, (2) 

where the passive sol and gel stresses t^gQi/^ggi are determined by 
linear constitutive laws and /soi/gei force densities. The active 
stress generated in the bulk of the gel phase is expected to be 
isotropic and hence, given by a multiple of the unit tensor: 
^gel ~^a^' assume only small deformations and thus restrict 
ourselves to linear elastic theory. We consider the gel phase as a 
porous viscoelastic active material [30] that is able to exert 
contractile stresses by interaction of the myosin-motor system with 
the actin filaments. The filament orientation in the cortex of 
Physarum that mainly determines the active and passive properties 



The factors y]^sof/gel front of the symmetric traceless parts 
denote the shear viscosities and yflfl^g^i the bulk viscosities for sol 
and gel phase, respectively. In the elastic stress G is the shear 
modulus and K the compression modulus. The mechanical force 
balance equations in the final form (for details, see Text S 1 in File 
S2) together with the incompressibility condition read: 

+ </ + iTa -p)\) + {ipl,(v -u) =0 
'^i&^:,-pl)-{ipl,{v-u) =0 (6) 

V-(P°,« + pLv) =0. 
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Model: Chemical part 

Calcium ions play a key role in the regulation of contraction in 
Physarum polycephalum. To describe the calcium kinetics, we use the 
biochemically realistic model presented by Smith and Saldana 
[31]. This model describes an oscillation mechanism that is driven 
by a phosphorylation-dephosphorylation cycle of myosin light 
chain kinases (MLCKs). A crucial feature of the Smith- Saldana 
model are autonomous calcium oscillations that occur already in 
the absence of mechanical feedback. Another entirely different 
mechanism was introduced to explain oscillations of plasmodial 
strands of Physarum, where the necessary feedback for the 
oscillations is provided by mechano-sensitive channels. This 
feedback acts upon a non-oscillatory calcium reaction kinetics 
and is therefore required to obtain the oscillations [12,39]. This 
model, however, predicts that the free calcium concentration and 
the mechanical tension oscillate in phase, whereas experiments 
exhibit an antiphase oscillation with a phase shift of n between 
these two quantities [33]. Moreover, experiments in the homog- 
enate of Physarum plasmodium wherein deformations and 
mechanical stresses are not possible yield calcium oscillations 
giving further evidence for an autonomous calcium oscillator [40] . 
The Smith-Saldana model [31] can be reduced to two ordinary 
differential equations involving the free calcium concentration ric 
and the fraction of phosphorylated myosin-light-chain kinases 

nc = {khiNc -rib- ric) kyric ^ - 0) + ^^0)^ / 

= :fc{nc4) 
j) = kQ(nc)(l-(t))-kE(l)= -./^{nc^)- 



ta = (FTeM)-Ta)/TT, 



(10) 



where Tt is the relaxation time the tension needs to approach its 
equilibrium value and Ft represents the mechanochemical 
coupling strength. According to the model in Ref. [31] the 
fraction of activated myosin motors (available for binding to actin) 
is given by 

Oa/bM = (kp(\-q2a/h(nc)))/(kD + kp(l-q2a/b(nc))) (11) 

qia/bM = {Ka/bflcf /{l+Ka/bflcf 



In this equation two additional parameters, the phosphorylation 
rate kp and dephosphorylation rate kj) of the LC 1 binding site 
come into play. The relaxation time constant Tj^ is introduced to 
obtain a phase shift between calcium concentration and active 
tension that is %. Earlier, a direct proportionality was used [31]. 
But this work could not reproduce the phase shift between calcium 
oscillation and active stress that was observed in experiments. 
Thus, we have adjusted Tt to obtain the correct phase relation. 

The next step in the derivation of the model is a spatial 
extension of Eq. (7) to a reaction-diffusion-advection system. The 
only species in this model that is transported by diffusion is the free 
calcium ric- Altogether, the following equation is obtained 



(12) 



These equations involve the myosin-bound calcium concentra- 
tion 



where w is the fluid velocity in the body-reference frame 
introduced above and Dc denotes the diffusion constant of the 
free calcium in the cytosol. 

Summary of the model 

Collecting all pieces of the model described so far, we obtain the 
following two-dimensional system of partial differential equations: 



n,M) = INm [ T^^il - ^) + T^^) 



and the calcium-dependent phosphorylation rate of the MLCK 



(8) 



+ ViTa-p) + pl,li{v-u) 



(13) 



1 r \ , u / K*nc 



(9) 



0 = »7*-Av + -jif V(V- V) - V;, - p» ;6(v - «) 



(14) 



The chemical parameters introduced in Ref [31] include the 
concentration of myosin Nm, the equilibrium calcium concentra- 
tion Nc, a pumping and leaking rate ky and kL of calcium 
vacuoles, and the dephosphorylation rate of the MLCK ks- 
Furthermore, we have the affmity Ka^b of calcium to the LCI 
binding site on the myosin light chain (MLC). The value of this 
affinity depends on whether the LC2 binding site is phosphory- 
lated (b) or not (a). The AC-cAMP-PKA signal cascade is modeled 
by Eq. (9) using a Hill-equation with coefficients ^, K^, and k^Q. 
For a detailed explanation of the terms in Eq. (7), we refer the 
reader to the publications [31,41]. The model is completed by a 
third equation that relates the activated fraction of myosin to the 
active tension Ta. In contrast to the original work of Saldana and 
Smith, we introduce a relaxation equation analogous to models of 
cardiac myocytes [42]: 



0 = V-(p»> + p«„,v) 



d,Ta = (FTe(nc,(l>)-Ta)lxT 



d,nc = - V-(nc{y - it)) + DcAric +fc{nc4) 



St<l>=Mnc,<l>)- 



(15) 



(16) 



(17) 



(18) 



Here, we have introduced the shear and bulk viscosities yi^sof/gel 
and y]^sofjggi of sol and gel phase and the linear elastic shear and 
compression modulus of the gel phase K and G, respectively. Note 
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Table 1. Mechanical parameters. 





Parameter 


Value 


Description 


Ft 


0-350 


mechanochemical coupling strength [45] 


Dc 


0.03 — 

min 


free calcium diffusion coefficient [43] 


K 


8.9 kPa 


gel compression modulus [45,46] 


G 


8.9 kPa 


gel shear modulus [45,46] 


shear 
^sol 


l.QPas 


effective sol shear viscosity [49,50] 


^shear 


l.OPas 


effective gel shear viscosity; indirectly from [51] 


^bulk 
^sol 


0.0 Pas 


effective sol bulk viscosity 


„bulk 
^gel 


0.0 Pas 


effective gel bulk viscosity 




5-10^ -5-10^ kg /{mm' min) 


drag coefficient, related to pore size [34,35] 




0.75 


sol volume fraction [44] 



Standard parameter set used for the mechanical part of the model given by Eqs. (13) - (15). 
doi:1 0.1 371/journal.pone.0099220.t001 



that the sol fraction and the free calcium concentration are given 
in the body-reference frame. The pressure p formally appears as a 
Lagrangian multiplier, due to the incompressibility condition (like 
in the incompressible Navier-Stokes equation). 

The above equations are defined in a circular domain that 
mimics the geometry of the Physarum droplets in experiments 
[18,19] (and U. Strachauer & M.J.B. Hauser, unpublished data). 
Zero displacement and zero cytosolic flow are imposed at the 
boundaries. Assuming, a non-permeable membrane of the droplet, 
no-flux boundary conditions \nc{x,t)'n{x,t) = () are employed for 
the calcium concentration Uc, where n is the normal vector at the 
boundary. As initial conditions we apply a small random noise ^ 
with zero mean as a perturbation to the homogeneous steady state 
solution. 



To compare with the experimentally measured height profiles of 
the droplet, one needs to estimate the local height field H{x,t) that 
does not appear explicitly in the two-dimensional model. We relate 
the relative height deviation h{x,t) : =(^(x,0 — ^o)/^o (where 
is the height in the reference configuration) to the divergence of 
the displacement field computed in the two-dimensional model: 



/z(jc,OocV-w(jc,0. 



(19) 



This approximation is based on the idea that deformations are 
locally isotropic [35]. Though Eqs. (13)-(18) may appear quite 
complex at first glance, their essential aspect is the combination of 



Table 2. Chemical parameters. 



Parameter 


Value 


Description 


kL 


0.24 min~^ 


leaking rate of vacuoles [31] 


ky 


4.8 min~^ 


pumping rate of vacuoles [31] 


k^ 
Q 


60 min~^ 


max. phosphorylation rate of MLCK, appears in 
function kqinc), [31] 


kE 


6.0 min~^ 


dephosphorylation rate of MLCK [31] 


kp 


30 min~^ 


phosphorylation rate of LC\ 


in function 9inc,(t>), [31] 


ko 


12 min ^ 


dephosphorylation rate of LC\ 


in function 9{nc,(t>), [31] 




1.5 nM-^ 


effective activation constant for the 


AC-cAMP-PKA chain, in function A:g(/2,.), [31] 


Ka 


2.3 nM-^ 


C<a!^"^ affinity with dephosphorylated MLCK 


Kb 


0.15 ^M-^ 


Cfl^+ affinity with phosphorylated MLCK [31] 


Nc 


25 /iM 


equilibrium total calcium concentration [31] 


Nm 


10 fiM 


total myosin concentration [31] 


It 


0.2 min 


relaxation time for tension generation [41] 



Standard parameter set used for the chemical part of the model given by Eqs. (16)-(18). 
doi:1 0.1 371/journal.pone.0099220.t002 
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the mechanical dynamics of a poroelastic medium and the calcium 
oscillator. 

Parameters 

A summary of the parameters used in the model is given in 
Tables 1 and 2. The values for the chemical oscillator are taken 
from Ref [31]. For the affinity Ka that appears in the functions fc 
and 6 in Eq. (7) two different values are considered. This 
parameter represents the affinity of the myosin light chain kinase 
to calcium ions in the unphosphorylated state. Autonomous 
calcium oscillations are obtained for the parameter Ka = 2.3 
liM~^, whereas a stationary calcium concentration is found for 
Ka = 2.0 iiM~^ in the absence of mechanical feedback. The effect 
of the parameter Ka on the autonomous calcium oscillator model 
for Physarum was studied systematically in Refs. [31,41]. For the 
diffusion coefficient of the free calcium we use a typical value of 
Dc = 0.03 mm^/min for small ions in cytoplasm [43]. 

In Physarum, the percentage of actin in the cytoplasm is 
estimated to be in the range 1 5 — 25% [44] . Since other proteins 
also contribute to the solid gel phase, we set p^^^ = 0.25 and 

p^^/ = 0.75. The Young modulus of a Physarum strand was 
determined to be about 10 kPa [45]. For sponge-like materials, 
measurements show that the Poisson ratio is very low: v^O [46] 
implying G^K in two dimensions. Therefore, we have set the 
parameters G = K=S.9 kPa. 

A typical value for the generated tension in plasmodial strands is 
Ta^20 kPa[47]. According to Eq. (16), the active tension Ta 
relaxes to an equilibrium value of T^ = FtO. Because 
^max<0-l(see Ref [48]), values up to ^^ = 350 kPa for Ft have 
been considered. Moreover, we have varied Ft in our study to 



illustrate the influence of the strength of mechanochemical 
coupling on the pattern dynamics. 

For the dynamic sol viscosity ^^^/^^ in Physarum, values in the 
range of 0.1— 0.5 Pas where measured [49]. Alternatively, one 
can calculate the sol viscosity from velocity profiles measured in 
Physarum [50]. With this method we obtain a dynamic sol 
viscosity of around 10 Pas. In the model, the value of the sol 
viscosity is then set to ^l^T^ ~ ^ Pas. In a study of a composite 
network containing actin filaments and microtubuli [51] the 
frequency dependent complex dynamic shear modulus 
G((jo) = G'(co) + iG"{w) was measured. For a viscoelastic material 
described by the Kelvin-Voigt model, one identifies 
G"(co) = (jorj^gfl^\ With a typical frequency (jo = 2nmin~^ of 
oscillations in Physarum a value of rj^^Jj'^^ ^ \Pas is obtained. For 
this value of gel viscosity the timescale, on which viscous gel stress 
is relevant is of the order of 10~^^. Assuming that the relevant 
timescale here is about Imin, viscous stress in the gel can be 
neglected. However, due to a stabilizing effect on the numerics, we 
keep this term. We neglect the influence from the bulk viscosities 
and set ^J"//gg/ = 0. Since the permeability kcci^^^^, the drag 
coefficient is P^^sol/^pore^ where ipore is the average pore size. 
Porous structures in Physarum exhibit several spatial scales: The 
actin network pore size is ipore = ^-25 fim [34] and the membrane 
structure of invaginations possesses pores with a typical diameter 
of about 10 fim [35]. If one assumes that larger porous structures 
determine the permeability, a drag coefficient of jS^^lO^ 
kg / (mrrr' min) is obtained. Since, however, the porous structure 
of the Physarum cytoskeleton may vary substantially over time, we 




Figure 3. Linear dispersion relation for varying mechanical coupling strength. Real part (left) and imaginary part (right) of the branch of 
the dispersion relation with largest real parts of the eigenvalues for a) a stable HSS {Ka = 2.0 nM~^) and b) an unstable HSS {Ka = 2.3 nM~^). The 
imaginary part is nonzero for all unstable wavenumbers: lm()i(q))^0. The mechanochemical coupling strength Ft is varied in each case a) and b) 
increasing from dark to light blue. The drag coefficient is chosen to be P = 5-10'^ kg/(mm^min) and the remaining parameters are given in Tables 1 
and 2. 

doi:1 0.1 371 /journal. pone.0099220.g003 



PLOS ONE I www.plosone.org 



6 



June 2014 | Volume 9 | Issue 6 | e99220 



A Model for Pattern Formation in Physarum Droplets 




< 



^ 2 



2 3 4 5 
wavenumber q/mm"^ 




o 




drag coefficient yS/(kg/(mm^min)) 

Figure 5. Phase diagram. Phase diagram in the plane spanned by 
the mechanochemical coupling strength Ft and the drag coefficient p. 
The black line denotes the threshold coupling strength F^^'', where the 
most unstable mode according to linear stability analysis of the HSS is 
nonzero. The dashed blue curve separates the plane according to Eq. 20 
in regions with Pe>l(larger Ft) and Pe<l(smaller Ft). 
doi:10.1371/journaLpone.0099220.g005 
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Figure 4. Linear dispersion relation for varying drag coeffi- 
cient. Branch of the dispersion relation with largest real part a) for 
different drag coefficients p that range over three orders of magnitude 
(logarithmic scale, from dark blue to cyan). The mechanochemical 
coupling strength is kept constant at Ft^2S kPa. b) Wavelength 
Ac = 2n/qc of the fastest growing mode versus drag coefficient p for 
three different mechanochemical coupling strengths Ft = 14 kPa (solid 
line), Ft = 44 kPa (dashed line) and Fr = 140 kPa (dotted line). The 
remaining parameters are given in Tables 1 and 2. 
doi:1 0.1 371 /journal. pone.0099220.g004 

have varied the parameter P in the range from 6.0-10^ to 

9.6.10^—^. 

mm^min 

Linear stability analysis 

Linear stability analysis is used here to investigate the 
spatiotemporal instability near the homogeneous steady state 
(HSS) without deformation and without any fluid motion. This 
HSS is given by = v^ = 0, = const, and the steady state values 
of chemical components and active tension that follow from the 
implicit equations /^(w^, = O,/^(n^,0'*') = 0 and = FTO{rf^,(j)^). 
The growth rate and the frequency (for imaginary eigenvalues) of a 
small perturbation [dfic, STa, Su, Sv, dp)e^^^^^^ to the HSS is 
given by the dispersion relation }^{q). 

Numerical integration 

The integration of the full PDE model given by Eqs. (13) - (18) 
was done with a hybrid method consisting of a finite element 
(FEM) and finite volume discretization scheme (FVM). This part 
was implemented as a stand-alone software written in C. Two- 



dimensional meshes, however, were generated with the free 
software Triangle [52]. With the operator-splitting technique the 
time integration step was divided into substeps that are carried out 
with different types of solvers: a fourth order Runge-Kutta method 
for the nonlinear reaction part, a linear FEM with implicit Euler 
stepping for the parabolic and elliptic parts and a FVM using the 
dual mesh of the triangulation (Voronoi diagram) for the advection 
step. (Operator splitting leads to subproblems of the form 
dtnc= Dl^nc-\- ...i^^vdiboYic PDE), a linear elastic problem 
(/xA + (;U + A) W-)u= ...(elliptic PDE) and an advection problem 
dtnc-\-^'{ncW) = ... (hyperbolic PDE).) The number of nodes for 
the mesh discretizing a disc with radius r = 1 mm was TV = 5218. 
This results in a typical mesh resolution of Ax 0.024 mm, 
whereas wavelengths obtained in the simulations do not go below 
^mw ^0.6 mm. The ratio Axj\min < 1/40 is assumed to provide a 
sufficient accuracy in order to resolve the profiles of the wave 
patterns and allow us to qualitatively compare our simulation 
results with experimental data. The time step size A/ was chosen to 
be 0.01 mz>2 that is about one order of magnitude smaller than the 
fastest time scale in the reaction kinetics given by Eqs. (16) - (18). 
The typical simulation length was 100 min. 

Results 

Linear stability analysis and dispersion relation 

Here, the mechanochemical coupling strength Ft and the drag 
coefficient P are varied to reveal their influence on the stability of 
the HSS and the shape of the dispersion curve A(g = |^|) . We have 
carried out numerical simulations where other model parameters 
were varied (results not shown) and found that indeed the drag 
coefficient P and the coupling strength Fj lead to the largest 
qualitative changes. In contrast, changes in sol viscosity by one 
order of magnitude had almost no effect on the simulation results. 
In Fig. 3 the branch of eigenvalues with the largest real part is 
shown for two scenarios (a and b) with different affinity Ka and for 
different values of Fj. The real and imaginary parts of the 
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Figure 6. Traveling wave, a) Snapshot of the free calcium 
concentration ric and c) relative height field h in color and the 
protoplasmic flow field v shown by arrows with length x |v|. Space-time 
plot of ric b) and h d) along the clotted line in panel a). The period of 
local oscillations is r=1.8 min. The parameters are Fr = 18 kPa and 
jS = 5.0-10^ kg /(mm^ min). The remaining values can be found in Tables 1 
and 2. A video file displaying the spatiotemporal dynamics including 
the transient initial phase corresponding to subfigure c) is included in 
the supporting material (see Movie SI in File SI). 
doi:1 0.1 371 /journal. pone.0099220.g006 

dispersion relation are displayed in separate plots. Fig. 3 a) shows 
the case where the HSS is stable in absence of mechanochemical 
coupling {Ft = 0). If Ft is increased above a critical value, the 
HSS exhibits a wave instability (oscillatory Turing instability). In 
Fig. 3 b) we consider a case where the HSS is already unstable 
against oscillations for Ft = 0, i.e. Re(k(0))>0. For larger 
wavenumbers q of the perturbation the imaginary part of the 
growth rate increases indicating that the frequency of waves is 
larger than the frequency of homogeneous oscillations. 

In the case shown in Fig. 3 b), there are two mechanisms that 
destabilize the HSS: the homogeneous oscillatory mode with q = 0 
that originates from the calcium kinetics described by Eqs. (7) and 
a wave-like perturbation at finite wavenumber (^t^O) induced by 
the mechanical feedback. Above a critical mechanochemical 
coupling strength Ft > F^^ > 0 the fastest growing mode (mode 
with largest real part of the eigenvalue) has a finite wavelength 
^7^0 and a nonzero imaginary part indicating wave dynamics. 

In Fig. 4 a) the influence of variation of the drag coefiicient P on 
the dispersion relation for fixed Ft is shown. The wavelength 
Kc = 2n/qc of the perturbation with largest growth rate increases 
with growing until the maximum with finite wavelength in the 
dispersion curve disappears. This discontinuity is visible in Fig. 4 b) 
where is plotted versus the drag coefficient ^ for different values 
of Ft. This figure also shows that the fastest growing wavelength 




Figure 7. Standing wave, a) Snapshot of the free calcium 
concentration ric and c) relative height field h in color and the 
protoplasmic flow field v shown by arrows with length x |v|. Space-time 
plot of ric b) and h d) along the clotted line in panel a). The period of 
local oscillations is r = 2.0 min. The parameters are Ft = 22 kPa and 
^ = 5.0-10'^ kg /(mm^ min). The remaining values can be found in Tables 1 
and 2. A video file corresponding to subfigure c) is included in the 
supporting material (see Movie S2 in File SI). 
doi:10.1371/journal.pone.0099220.g007 

Ac decreases with the mechanochemical coupling strength. 
Nevertheless, the wavelength depends only weakly on over a 
large range of values. 

The linear stability analysis gives already an important insight 
into the essential physics of the model defined by Eqs. (13) - (18): 
The mechanical part leads to a fastest growing mode with non- 
zero wavenumber while the calcium oscillator can switch the long- 
wavelength instability observed in the purely mechanical system 
[30] to a short wavelength instability by removing the conservation 
constraint for the integrated concentration of free calcium. 

Numerical simulations 

For different values of the coupling strength Ft and the drag 
coefficient fi, we present numerical simulations of the dynamics of 
full nonlinear model equations. The corresponding phase diagram 
is shown in Fig. 5. Therein, the solid black line separates the phase 
plane according to the shape of the dispersion curve obtained from 
the linear stability analysis of the HSS: It shows the threshold value 
as a function of the parameter p. When Ft < Fj^ there is no 
discrete wavenumber qk>(^, for which X(q/^)>X(qQ=0). In the 
opposite case F > Fj^ there exists at least one qk>^ with 
X(qk)>X(qo). The prediction of homogeneous oscillations for 
Ft < Fj^ from linear stability analysis is confirmed by numerical 
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Figure 8. Single rotating spiral, a) Snapshot of the free calcium 
concentration ric and c) relative height field h in color and the 
protoplasmic flow field v shown by arrows with length x |v|. Space-time 
plot of He b) and h d) along a circle marked by the dotted line in a). The 
period of local oscillations is T^l.^min. The parameters are Ft = 22 
kPa and P = 5.0-10'^ kg/(mtrr'min). For the remaining values see Tables 1 
and 2. A video file corresponding to subfigure c) is included in the 
supporting material (see Movie S3 in File SI). 
doi:1 0.1 371 /journal. pone.0099220.g008 

simulations for P<10^ kg / {mn? min) (see Fig. 5, blue discs). 
However, for larger values of P one finds a considerable 
discrepancy between linear stability analysis and simulation results. 
For homogeneous oscillations in the calcium concentration the 
flow and deformation fields both vanish. 

To address the question how much influence the advective 
coupling relative to diffusion has, one can consider a Peclet 
number that is given by the ratio of diffusive to advective time 
scales (see also Ref [22] for a motivation of this definition) 



Pe = e,„a.FTl{D,P), 



(20) 



where dmax^^-^^ is the typical amplitude for the oscillations in 
the variable 0 appearing in Eq. (10) for Ka = 23 fiM'K In Fig. 5 
the blue line corresponds to Pe=\. Note, that above the line 
Pe>l, there are no homogeneous oscillations and different types 
of patterns prevail. In some cases, simple patterns like rotating 
spirals are traveling waves are also found for Pe<l. Altogether, 
this consideration shows that the mechanochemical coupling has 
to be strong enough to overcome the homogenizing effect of 
diffusion for mechanochemical waves and patterns to emerge. 

Above F^r^!^, traveling, standing and spiral waves occur in the 
simulations. In Fig. 6 a traveling wave is depicted by snapshots and 
space-time plots. In addition to the free calcium concentration ric 




y-:.^^ v.^^^ C':-- 



Figure 9. Radial wave, a) Snapshot of the free calcium concentration 
Yic and c) relative height field h in color and the protoplasmic flow field v 
shown by arrows with length x|v|. Space-time plot of Uc b) and h d) 
along the dotted line in a). The period of local oscillations is T=lJmin. 
The parameters are = 194 kPa and ^ = 5.0-10^ kg/(mm^min). For the 
remaining values see Tables 1 and 2. A video file corresponding to 
subfigure c) is included in the supporting material (see Movie S4 in File 
SI). 

doi:10.1371/journal.pone.0099220.g009 

the plots show the relative height field h and the sol velocity v, 
since these quantities are most likely to be measured in 
experiments. We get an average velocity of the wave front that 
is about 0.086 mm/s. 

For values of ^ around 5-10^ kg / (mnr' min) standing waves are 
observed for the relative height field h, while for the concentration 
field He a traveling wave appears (see Fig. 7). 

Traveling and standing waves often coexist with single rotating 
spiral patterns. The pattern selection depends on the initial 
conditions, in particular on the number of phase singularities. Two 
counter-rotating spirals annihilate and give way to a traveling 
wave, whereas a single rotating spiral is stable. A single rotating 
spiral is presented in Fig. 8. No distinct direction of propagation or 
rotation is preferred reflecting the symmetries of the system. The 
coexistence region of traveling waves and spirals is marked by 
violet squares in the phase diagram in Fig. 5. From the space-time 
plots in Fig. 8 c) one finds a wave speed of 0.047 mm/s. 

For larger mechanochemical coupling strength Fj^ there is no 
coexistence of spirals with traveling waves: a rotating wave is the 
only attractor (see Fig. 5, green triangles). For large enough values 
of the drag coefficient ^ a mode with largest possible wavenumber 
Ai = 2L determines the emerging patterns. This is confirmed by 
the numerical simulations. For jS<10^ kg/(mm^min) periodic 
patterns are obtained, even for large coupling Ft (see Fig. 5, 
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0.34 — ^ space/mm 




Figure 10. Irregular wave pattern, a) Snapshot of the free calcium 
concentration ric and c) relative height field h in color and the 
protoplasmic flow field v shown by arrows with length x |v|. Space-time 
plot of He b) and h d) along the dotted line in a). The parameters are 
Ft ^356 kPa and ^ = 5.0-10^ kg/(mm^min). For remaining values see 
Tables 1 and 2. A video file corresponding to subfigure c) is included in 
the supporting material (see Movie S5 in File SI). 
doi:10.1371/journal.pone.0099220.g010 

brown squares). These patterns have a characteristic wavelength of 
the same order as the system size. The simplest pattern is an 
antiphase oscillations in the form of a radial wave that is reflected 
at the boundaries (see Fig. 9) that has the symmetry as the 
experimentally found pattern in Fig. Id. 

For large drag coefiicients (see Fig. 10 and phase diagram in 
Fig. 5, pink circles), irregular patterns like the experimental pattern 
in Fig. lb with a wavelength significantly smaller than the system 
size are obtained. For the wave segments in these spatiotemporal 
patterns we get a typical velocity of 0.03 mm/ s that is much slower 
than that of traveling or spiral wave. 

Discussion 

In this article we have combined a recently published novel 
mechanical continuum model of the cytoplasm as a poroelastic 
active medium [30] with the Smith-Saldana model for calcium 
oscillations in Physarum protoplasm [31]. Upon increase of the 
mechanochemical coupling strength the homogeneous steady state 
in this model is destabilized by an oscillatory Turing instability 
with finite wave number. Homogeneous oscillations of the calcium 
concentration are replaced by spatiotemporal patterns and waves 
connected with local deformation and fluid motion. Practically all 



spatiotemporal contraction patterns observed experimentally in 
Refs. [18,19] including rotating spirals, traveling and standing 
waves, antiphase oscillation and irregular, chaotic waves could be 
reproduced by numerical simulations of this model. In contrast to 
earlier models (see e.g. Ref [20]) for Physarum protoplasm, a 
closed set of mechanical force-balance equations is derived that 
allows for explicit computations of pressure and flow fields. The 
cytoplasm is treated as a two-phase material, consisting of a passive 
fluid sol phase and an active solid viscoelastic gel phase. The basic 
idea was already sketched in [41], where however only the case 
was considered, where the mechanochemical coupling is approx- 
imated by a global coupling in the dynamics for calcium. Here, we 
have instead analyzed and simulated the interplay of mechanical 
deformation of the cytoskeleton, fluid flow in the cytosol and 
chemical (= calcium) concentration and the associated spatio- 
temporal dynamics. 

Unlike other models that treat the cytoskeleton as a viscoelastic 
fluid (see e.g. Ref [53]), we consider the cytoskeleton to be a 
viscoelastic solid described the Kelvin-Voigt model. This is a valid 
approximation in Physarum, since the relaxation time of elastic 
tension is larger by about a factor of three than the calcium 
oscillation period in the system [32]. The porosity of the solid gel 
phase allows a flow of cytosol that carries calcium and regulates 
the tension generation. This is represented by an advective 
transport of calcium in addition to diffusion and introduces a 
nonlinear feedback mechanism that couples flows and deforma- 
tions to a nonlinear reaction kinetics. Related models of active gels 
and fluids [22,24,54] consider the transport of motors, whereas we 
have considered the transport of calcium that acts as a motor- 
regulating species in Physarum. 

We have limited our considerations to small displacements u 
and deformations |Vi/|. Thus, the equations are formulated up to 
linear order in the displacement field and its gradient. 

A linear stability analysis of a homogeneous steady state with 
constant calcium concentration and zero deformation and flow has 
been carried out for the presented model. The resulting dispersion 
relations reveal that the mechanical feedback provides a new 
mechanism of an oscillatory Turing instability with nonzero 
wavenumber if the mechanochemical coupling strength is positive 
Ft>0. In numerical simulations, homogeneous oscillations get 
destabilized for sufficiently large coupling strength Ft- Rotating 
spirals and traveling waves obtained in the simulations are 
common patterns found in experiments [18,19] (and U. 
Strachauer & M.J.B. Hauser, unpublished data). The wave speed 
we obtain for traveling and spiral waves agree with the findings in 
[18]. Note, however, that the antiphase and irregular patterns 
(Figs. 9 and 10) were obtained at large coupling strength Ft 
yielding also large deformations. Thus, these results violate the 
small deformation assumption of our model and should be used 
only for qualitative comparisons. For a quantitative treatment of 
large deformations, an extension of the above model to a nonlinear 
elasticity formulation [55] becomes necessary. 

The mechanism of the experimentally observed transition 
between different patterns on a timescale much larger than the 
typical calcium oscillation period is not fully understood. The 
phase diagram in Fig. 5 shows that a minor change in the coupling 
strength Ft or drag coefffcient P can result in a different type of 
pattern. Hence, a variation of one of these parameters over the 
time of the experiment may be responsible for qualitative changes. 

For the sake of simplicity, we choose a simple Kelvin-Voigt 
model as a starting point and found that it is sufficient to 
reproduce the experimental findings on deformation patterns in 
Physarum. Future work, e. g. on the modeling of moving droplets 
of Physarum, may require a more complicated model that exhibits 
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solid behavior of the cytoskeleton at short time scales and fluid 
behavior at long time scales. 

In order to model the stage of migration of a Physarum droplet 
free-boundary conditions have to be imposed and an interaction 
with the substrate must be considered (see e.g. [36]). Such an 
extension of the model may eventually help to understand the 
interplay of chemical and mechanical processes in the self- 
organized amoeboid movement of Physarum droplets. 

Supporting Information 

Comments about the sol continuity equation, the osmotic 
pressure, as well as a derivation of the force balance equations 
from an extremal principle can be found in the Supporting 
Information (Text SI in File S2) provided with this article. 
Furthermore, we include movies (see Movies S1-S5 in File SI) with 
the spatiotemporal dynamics corresponding to Figs. 6-10. 

Supporting Information 

File SI Supporting movies. Movie SI. Spatiotemporal 
dynamics of traveling wave pattern: Relative height field h in 
color and flow field v given by arrows corresponding to Fig. 6. The 
simulation time is 100mm beginning with the initial state. Movie 
S2. Spatiotemporal dynamics of standing wave pattern: Relative 
height field h in color and flow field v given by arrows 
corresponding to Fig. 7. The simulation time is 100 min beginning 
with the initial state. Movie S3. Spatiotemporal dynamics of 
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